R Notebook
library(ggplot2)
library(ggsci)
library(hdf5r)
library(patchwork)
library(RColorBrewer)
library(Seurat)
library(tidyverse)setwd("/Users/christinacomo/Desktop/Spatial")
getwd()[1] "/Users/christinacomo/Desktop/Spatial"data_dir <- '/Users/christinacomo/Desktop/Spatial/data/B'
list.files(data_dir) # Should show filtered_feature_bc_matrix.h5[1] "cellshighlighted_b.png" "Clusterandspatial_B.png" "filtered_feature_bc_matrix.h5"
[4] "nCount_B.png" "spatial" "SpatialHighlited_B.png"
[7] "web_summary.html" # grab the spot barcodes to use for subsetting
spots.use <- regions_of_interest$barcodes
regions_of_interest <- regions_of_interest %>% column_to_rownames(var = "barcodes")sample <- "control2"
SampleB[["orig.ident"]] <- sample
SampleB <- SetIdent(SampleB, value = so@meta.data$orig.ident)SampleB <- PercentageFeatureSet(SampleB, "^mt-", col.name = "percent_mito")
SampleB <- PercentageFeatureSet(SampleB, "^Hb.*-", col.name = "percent_hb")VlnPlot(SampleB,
features = c("nCount_Spatial",
"nFeature_Spatial",
"percent_mito",
"percent_hb"),
pt.size = 0.1,
ncol = 2) + NoLegend()
SpatialFeaturePlot(SampleB,
features = c("nCount_Spatial",
"nFeature_Spatial",
"percent_mito",
"percent_hb"))
SampleB <- SCTransform(SampleB, assay = "Spatial")Calculating cell attributes from input UMI matrix: log_umi
Variance stabilizing transformation of count matrix of size 16658 by 1620
Model formula is y ~ log_umi
Get Negative Binomial regression parameters per gene
Using 2000 genes, 1620 cells
|
| | 0%
|
|=========================== | 25%
|
|====================================================== | 50%
|
|================================================================================ | 75%
|
|===========================================================================================================| 100%There are 3 estimated thetas smaller than 1e-07 - will be set to 1e-07
Found 108 outliers - those will be ignored in fitting/regularization step
Second step: Get residuals using fitted parameters for 16658 genes
|
| | 0%
|
|=== | 3%
|
|====== | 6%
|
|========= | 9%
|
|============= | 12%
|
|================ | 15%
|
|=================== | 18%
|
|====================== | 21%
|
|========================= | 24%
|
|============================ | 26%
|
|=============================== | 29%
|
|=================================== | 32%
|
|====================================== | 35%
|
|========================================= | 38%
|
|============================================ | 41%
|
|=============================================== | 44%
|
|================================================== | 47%
|
|====================================================== | 50%
|
|========================================================= | 53%
|
|============================================================ | 56%
|
|=============================================================== | 59%
|
|================================================================== | 62%
|
|===================================================================== | 65%
|
|======================================================================== | 68%
|
|============================================================================ | 71%
|
|=============================================================================== | 74%
|
|================================================================================== | 76%
|
|===================================================================================== | 79%
|
|======================================================================================== | 82%
|
|=========================================================================================== | 85%
|
|============================================================================================== | 88%
|
|================================================================================================== | 91%
|
|===================================================================================================== | 94%
|
|======================================================================================================== | 97%
|
|===========================================================================================================| 100%Computing corrected count matrix for 16658 genes
|
| | 0%
|
|=== | 3%
|
|====== | 6%
|
|========= | 9%
|
|============= | 12%
|
|================ | 15%
|
|=================== | 18%
|
|====================== | 21%
|
|========================= | 24%
|
|============================ | 26%
|
|=============================== | 29%
|
|=================================== | 32%
|
|====================================== | 35%
|
|========================================= | 38%
|
|============================================ | 41%
|
|=============================================== | 44%
|
|================================================== | 47%
|
|====================================================== | 50%
|
|========================================================= | 53%
|
|============================================================ | 56%
|
|=============================================================== | 59%
|
|================================================================== | 62%
|
|===================================================================== | 65%
|
|======================================================================== | 68%
|
|============================================================================ | 71%
|
|=============================================================================== | 74%
|
|================================================================================== | 76%
|
|===================================================================================== | 79%
|
|======================================================================================== | 82%
|
|=========================================================================================== | 85%
|
|============================================================================================== | 88%
|
|================================================================================================== | 91%
|
|===================================================================================================== | 94%
|
|======================================================================================================== | 97%
|
|===========================================================================================================| 100%Calculating gene attributes
Wall clock passed: Time difference of 18.19644 secs
Determine variable features
Place corrected count matrix in counts slot
Centering data matrix
|
| | 0%
|
|=========================== | 25%
|
|====================================================== | 50%
|
|================================================================================ | 75%
|
|===========================================================================================================| 100%
Set default assay to SCTSampleB <- RunPCA(SampleB, assay = "SCT", verbose = FALSE)
SampleB <- FindNeighbors(SampleB, reduction = "pca", dims = 1:30, verbose = FALSE)
SampleB <- FindClusters(SampleB, verbose = FALSE)
SampleB <- RunUMAP(SampleB, reduction = "pca", dims = 1:30)14:09:17 UMAP embedding parameters a = 0.9922 b = 1.112
14:09:17 Read 1620 rows and found 30 numeric columns
14:09:17 Using Annoy for neighbor search, n_neighbors = 30
14:09:17 Building Annoy index with metric = cosine, n_trees = 50
0% 10 20 30 40 50 60 70 80 90 100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
14:09:17 Writing NN index file to temp file /var/folders/2b/gcm1xy2n6s50jqf9f3ln87140000gn/T//RtmpWq9Udx/file58a42d3bd9cb
14:09:17 Searching Annoy index using 1 thread, search_k = 3000
14:09:17 Annoy recall = 100%
14:09:18 Commencing smooth kNN distance calibration using 1 thread
14:09:19 Initializing from normalized Laplacian + noise
14:09:19 Commencing optimization for 500 epochs, with 57888 positive edges
0% 10 20 30 40 50 60 70 80 90 100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
14:09:21 Optimization finishedp1 <- DimPlot(SampleB, reduction = "umap", label = TRUE)p2 <- SpatialDimPlot(SampleB, label = TRUE, label.size = 3)Scale for 'fill' is already present. Adding another scale for 'fill', which will replace the existing scale.p1 + p2
saveRDS(SampleB, "/Users/christinacomo/Desktop/Spatial/RDS_files/SampleB_initial.RDS")#Now subset with spots
rawdata_dir <- ('/Users/christinacomo/Desktop/Spatial/Raw Data')
MGE <- read.csv("/Users/christinacomo/Desktop/Spatial/Raw Data/Control2MGE.csv")
LGE <- read.csv("/Users/christinacomo/Desktop/Spatial/Raw Data/Control2LGE.csv")
Midline <- read.csv("/Users/christinacomo/Desktop/Spatial/Raw Data/Control2Midline.csv")
Dorsal_Cortex <- read.csv("/Users/christinacomo/Desktop/Spatial/Raw Data/Control2DorsalCortex.csv")
Ventral_Cortex <- read.csv("/Users/christinacomo/Desktop/Spatial/Raw Data/Control2VentralCortex.csv")# rename columns in cortex and meninges dataframes for merging
names(MGE) <- c("barcodes", "region")
names(LGE) <- c("barcodes", "region")
names(Midline) <- c("barcodes", "region")
names(Dorsal_Cortex) <- c("barcodes", "region")
names(Ventral_Cortex) <- c("barcodes", "region")# merge cortex and meninges dataframes together
regions_of_interest <- rbind(MGE, LGE, Midline, Dorsal_Cortex, Ventral_Cortex) # grab the spot barcodes to use for subsetting
spots.use <- regions_of_interest$barcodes
regions_of_interest <- regions_of_interest %>% column_to_rownames(var = "barcodes")# check dimensions of regions_of_interest and make sure it is the sum of cc_cortex and cc_meninges
dim(regions_of_interest)[1] 174 1# now, subset the seurat object
SampleB <- subset(SampleB, cells = spots.use)subset_B <- AddMetaData(SampleB, metadata = regions_of_interest)
p1 <- SpatialDimPlot(subset_B)p1
SpatialFeaturePlot(subset_B,
features = c("nCount_Spatial",
"nFeature_Spatial",
"percent_mito",
"percent_hb"))
saveRDS(subset_B, "/Users/christinacomo/Desktop/Spatial/RDS_files/SampleB_subsetwithbarcodes.RDS")LS0tCnRpdGxlOiAiUiBOb3RlYm9vayIKb3V0cHV0OiBodG1sX25vdGVib29rCmVkaXRvcl9vcHRpb25zOiAKICBjaHVua19vdXRwdXRfdHlwZTogaW5saW5lCi0tLQoKLS0tCnRpdGxlOiAiUiBOb3RlYm9vayIKb3V0cHV0OiBodG1sX25vdGVib29rCmVkaXRvcl9vcHRpb25zOiAKICBjaHVua19vdXRwdXRfdHlwZTogaW5saW5lCi0tLQpgYGB7cn0KbGlicmFyeShnZ3Bsb3QyKQpsaWJyYXJ5KGdnc2NpKQpsaWJyYXJ5KGhkZjVyKQpsaWJyYXJ5KHBhdGNod29yaykKbGlicmFyeShSQ29sb3JCcmV3ZXIpCmxpYnJhcnkoU2V1cmF0KQpsaWJyYXJ5KHRpZHl2ZXJzZSkKYGBgCgoKYGBge3Igc2V0dXAsIGluY2x1ZGU9RkFMU0V9CmtuaXRyOjpvcHRzX2NodW5rJHNldCh3YXJuaW5nID0gRkFMU0UsIG1lc3NhZ2UgPSBGQUxTRSwgZXJyb3IgPSBGQUxTRSkgCmBgYAoKCmBgYHtyfQpzZXR3ZCgiL1VzZXJzL2NocmlzdGluYWNvbW8vRGVza3RvcC9TcGF0aWFsIikKZ2V0d2QoKQpgYGAKCmBgYHtyfQpkYXRhX2RpciA8LSAnL1VzZXJzL2NocmlzdGluYWNvbW8vRGVza3RvcC9TcGF0aWFsL2RhdGEvQicKbGlzdC5maWxlcyhkYXRhX2RpcikgIyBTaG91bGQgc2hvdyBmaWx0ZXJlZF9mZWF0dXJlX2JjX21hdHJpeC5oNQpgYGAKCgpgYGB7cn0KU2FtcGxlQiA8LSBMb2FkMTBYX1NwYXRpYWwoCiAgZGF0YV9kaXIsCiAgZmlsZW5hbWUgPSAiZmlsdGVyZWRfZmVhdHVyZV9iY19tYXRyaXguaDUiLAogIGFzc2F5ID0gIlNwYXRpYWwiLAogIHNsaWNlID0gInNsaWNlMSIsCiAgZmlsdGVyLm1hdHJpeCA9IFRSVUUsCiAgdG8udXBwZXIgPSBGQUxTRSwKICBpbWFnZSA9IE5VTEwsCikKYGBgCgpgYGB7cn0Kc2FtcGxlIDwtICJjb250cm9sMiIKU2FtcGxlQltbIm9yaWcuaWRlbnQiXV0gPC0gc2FtcGxlClNhbXBsZUIgPC0gU2V0SWRlbnQoU2FtcGxlQiwgdmFsdWUgPSBzb0BtZXRhLmRhdGEkb3JpZy5pZGVudCkKYGBgCgoKYGBge3J9ClNhbXBsZUIgPC0gUGVyY2VudGFnZUZlYXR1cmVTZXQoU2FtcGxlQiwgIl5tdC0iLCBjb2wubmFtZSA9ICJwZXJjZW50X21pdG8iKQpTYW1wbGVCIDwtIFBlcmNlbnRhZ2VGZWF0dXJlU2V0KFNhbXBsZUIsICJeSGIuKi0iLCBjb2wubmFtZSA9ICJwZXJjZW50X2hiIikKYGBgCgpgYGB7ciwgZmlnLmhlaWdodCA9IDMsIGZpZy53aWR0aCA9IDN9ClZsblBsb3QoU2FtcGxlQiwgCiAgICAgICAgZmVhdHVyZXMgPSBjKCJuQ291bnRfU3BhdGlhbCIsIAogICAgICAgICAgICAgICAgICAgICAibkZlYXR1cmVfU3BhdGlhbCIsIAogICAgICAgICAgICAgICAgICAgICAicGVyY2VudF9taXRvIiwgCiAgICAgICAgICAgICAgICAgICAgICJwZXJjZW50X2hiIiksCiAgICAgICAgcHQuc2l6ZSA9IDAuMSwgCiAgICAgICAgbmNvbCA9IDIpICsgTm9MZWdlbmQoKQpgYGAKCmBgYHtyfQpTcGF0aWFsRmVhdHVyZVBsb3QoU2FtcGxlQiwgCiAgICAgICAgICAgICAgICAgICBmZWF0dXJlcyA9IGMoIm5Db3VudF9TcGF0aWFsIiwgCiAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgIm5GZWF0dXJlX1NwYXRpYWwiLCAKICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAicGVyY2VudF9taXRvIiwgCiAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgInBlcmNlbnRfaGIiKSkKYGBgCgpgYGB7cn0KU2FtcGxlQiA8LSBTQ1RyYW5zZm9ybShTYW1wbGVCLCBhc3NheSA9ICJTcGF0aWFsIikKCmBgYAoKYGBge3J9ClNhbXBsZUIgPC0gUnVuUENBKFNhbXBsZUIsIGFzc2F5ID0gIlNDVCIsIHZlcmJvc2UgPSBGQUxTRSkKU2FtcGxlQiA8LSBGaW5kTmVpZ2hib3JzKFNhbXBsZUIsIHJlZHVjdGlvbiA9ICJwY2EiLCBkaW1zID0gMTozMCwgdmVyYm9zZSA9IEZBTFNFKQpTYW1wbGVCIDwtIEZpbmRDbHVzdGVycyhTYW1wbGVCLCB2ZXJib3NlID0gRkFMU0UpClNhbXBsZUIgPC0gUnVuVU1BUChTYW1wbGVCLCByZWR1Y3Rpb24gPSAicGNhIiwgZGltcyA9IDE6MzApCmBgYAoKYGBge3J9CnAxIDwtIERpbVBsb3QoU2FtcGxlQiwgcmVkdWN0aW9uID0gInVtYXAiLCBsYWJlbCA9IFRSVUUpCnAyIDwtIFNwYXRpYWxEaW1QbG90KFNhbXBsZUIsIGxhYmVsID0gVFJVRSwgbGFiZWwuc2l6ZSA9IDMpCnAxICsgcDIKYGBgCgpgYGB7cn0Kc2F2ZVJEUyhTYW1wbGVCLCAiL1VzZXJzL2NocmlzdGluYWNvbW8vRGVza3RvcC9TcGF0aWFsL1JEU19maWxlcy9TYW1wbGVCX2luaXRpYWwuUkRTIikKYGBgCgoKYGBge3J9CiNOb3cgc3Vic2V0IHdpdGggIHNwb3RzIApyYXdkYXRhX2RpciA8LSAoJy9Vc2Vycy9jaHJpc3RpbmFjb21vL0Rlc2t0b3AvU3BhdGlhbC9SYXcgRGF0YScpCgpNR0UgPC0gcmVhZC5jc3YoIi9Vc2Vycy9jaHJpc3RpbmFjb21vL0Rlc2t0b3AvU3BhdGlhbC9SYXcgRGF0YS9Db250cm9sMk1HRS5jc3YiKQpMR0UgPC0gcmVhZC5jc3YoIi9Vc2Vycy9jaHJpc3RpbmFjb21vL0Rlc2t0b3AvU3BhdGlhbC9SYXcgRGF0YS9Db250cm9sMkxHRS5jc3YiKQpNaWRsaW5lIDwtIHJlYWQuY3N2KCIvVXNlcnMvY2hyaXN0aW5hY29tby9EZXNrdG9wL1NwYXRpYWwvUmF3IERhdGEvQ29udHJvbDJNaWRsaW5lLmNzdiIpCkRvcnNhbF9Db3J0ZXggPC0gcmVhZC5jc3YoIi9Vc2Vycy9jaHJpc3RpbmFjb21vL0Rlc2t0b3AvU3BhdGlhbC9SYXcgRGF0YS9Db250cm9sMkRvcnNhbENvcnRleC5jc3YiKQpWZW50cmFsX0NvcnRleCA8LSByZWFkLmNzdigiL1VzZXJzL2NocmlzdGluYWNvbW8vRGVza3RvcC9TcGF0aWFsL1JhdyBEYXRhL0NvbnRyb2wyVmVudHJhbENvcnRleC5jc3YiKQpgYGAKCmBgYHtyfQojIHJlbmFtZSBjb2x1bW5zIGluIGNvcnRleCBhbmQgbWVuaW5nZXMgZGF0YWZyYW1lcyBmb3IgbWVyZ2luZwpuYW1lcyhNR0UpIDwtIGMoImJhcmNvZGVzIiwgInJlZ2lvbiIpCm5hbWVzKExHRSkgPC0gYygiYmFyY29kZXMiLCAicmVnaW9uIikKbmFtZXMoTWlkbGluZSkgPC0gYygiYmFyY29kZXMiLCAicmVnaW9uIikKbmFtZXMoRG9yc2FsX0NvcnRleCkgPC0gYygiYmFyY29kZXMiLCAicmVnaW9uIikKbmFtZXMoVmVudHJhbF9Db3J0ZXgpIDwtIGMoImJhcmNvZGVzIiwgInJlZ2lvbiIpCmBgYAoKYGBge3J9CiMgbWVyZ2UgY29ydGV4IGFuZCBtZW5pbmdlcyBkYXRhZnJhbWVzIHRvZ2V0aGVyCnJlZ2lvbnNfb2ZfaW50ZXJlc3QgPC0gcmJpbmQoTUdFLCBMR0UsIE1pZGxpbmUsIERvcnNhbF9Db3J0ZXgsIFZlbnRyYWxfQ29ydGV4KSAKYGBgCgoKYGBge3J9CiMgZ3JhYiB0aGUgc3BvdCBiYXJjb2RlcyB0byB1c2UgZm9yIHN1YnNldHRpbmcKc3BvdHMudXNlIDwtIHJlZ2lvbnNfb2ZfaW50ZXJlc3QkYmFyY29kZXMKcmVnaW9uc19vZl9pbnRlcmVzdCA8LSByZWdpb25zX29mX2ludGVyZXN0ICU+JSBjb2x1bW5fdG9fcm93bmFtZXModmFyID0gImJhcmNvZGVzIikKYGBgCgpgYGB7cn0KIyBjaGVjayBkaW1lbnNpb25zIG9mIHJlZ2lvbnNfb2ZfaW50ZXJlc3QgYW5kIG1ha2Ugc3VyZSBpdCBpcyB0aGUgc3VtIG9mIGNjX2NvcnRleCBhbmQgY2NfbWVuaW5nZXMKZGltKHJlZ2lvbnNfb2ZfaW50ZXJlc3QpCmBgYAoKYGBge3J9CiMgbm93LCBzdWJzZXQgdGhlIHNldXJhdCBvYmplY3QKU2FtcGxlQiA8LSBzdWJzZXQoU2FtcGxlQiwgY2VsbHMgPSBzcG90cy51c2UpCmBgYAoKYGBge3J9CnN1YnNldF9CIDwtIEFkZE1ldGFEYXRhKFNhbXBsZUIsIG1ldGFkYXRhID0gcmVnaW9uc19vZl9pbnRlcmVzdCkKcDEgPC0gU3BhdGlhbERpbVBsb3Qoc3Vic2V0X0IpCnAxCmBgYAoKCmBgYHtyfQpTcGF0aWFsRmVhdHVyZVBsb3Qoc3Vic2V0X0IsIAogICAgICAgICAgICAgICAgICAgZmVhdHVyZXMgPSBjKCJuQ291bnRfU3BhdGlhbCIsIAogICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICJuRmVhdHVyZV9TcGF0aWFsIiwgCiAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgInBlcmNlbnRfbWl0byIsIAogICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICJwZXJjZW50X2hiIikpCmBgYAoKYGBge3J9CnNhdmVSRFMoc3Vic2V0X0IsICIvVXNlcnMvY2hyaXN0aW5hY29tby9EZXNrdG9wL1NwYXRpYWwvUkRTX2ZpbGVzL1NhbXBsZUJfc3Vic2V0d2l0aGJhcmNvZGVzLlJEUyIpCmBgYAoKCg==